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Abstract. Static and dynamical structure factors for shear and longitudinal modes 
of the velocity and density fields are computed for a granular system fluidized by a 
stochastic bath with friction. Analytical expressions are obtained through fluctuating 
hydrodynamics and are successfully compared with numerical simulations up to a 
C^ ' volume fraction ~ 50%. Hydrodynamic noise is the sum of external noise due to 

the bath and internal one due to collisions. Only the latter is assumed to satisfy the 
, fluctuation-dissipation relation with the average granular temperature. 

f-H . Static velocity structure factors S±{k) and S\\{k) display a general non-constant 

O ! behavior with two plateaux at large and small k, representing the granular temperature 

Tg and the bath temperature Tf, > Tg respectively. From this behavior, two different 
velocity correlation lengths are measured, both increasing as the packing fraction is 
raised. This growth of spatial order is in agreement with the behaviour of dynamical 
structure factors, the decay of which becomes slower and slower at increasing density. 
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1. Introduction 

Granular media display a wide catalog of non-equilibrium phenomena [1]. These 
materials are constituted by a number A^ of elementary constituents, grains of typical 
diameter between 0.1 and 10 mm. The number A^ ^ 1 is usually large enough to allow, 
or require, a statistical treatment. Unfortunately, interactions are non-conservative, 
resulting in the failure of equilibrium statistical mechanics. Kinetic theories, from 
Boltzmann equation to hydrodynamics [2], together with numerical simulations [3], are 
the best tools to describe those systems and to compare with real experiments, with the 
caveat of a proper adaptation to the peculiarity of granular interactions. 

One of the debated points of granular kinetic theories is the way noise should be 
added to hydrodynamics in order to describe mesoscopic fluctuations [U [5]. This is a 
general problem in non-equilibrium systems |6] (e.g. sheared fluids [7]), but here is even 
more pressing, given the rather small number of particles in a granular system: one 
has typically A^ ~ 10'^ -^ 10^, even in experiments, so that fluctuations can hardly be 
neglected. Moreover, in granular systems the dynamics is non- conservative and therefore 
Fluctuation-Dissipation relations do not hold in general [H |9l [TOl [TT], with exceptions 
in driven dilute cases [12], [131 E] • In non-dilute systems, it is therefore difficult even to 
define a temperature, making tricky the modelization of fiuctuations [151 [16] . 

A comprehensive study of the fiuctuating hydrodynamics of a driven granular fiuid 
is presented here. Static and dynamical structure factors are computed analytically 
in the framework of linearized hydrodynamics, and compared with extensive numerical 
simulations. A very good agreement is found between analytical and numerical results 
in a wide range of parameters, implying that, for this kind of model, fiuctuating 
hydrodynamics is able to describe large scale fiuctuations in a satisfactory manner. 

The peculiarity of the model we have studied, when compared to others present 
in the hterature [I?] [H], is the prescription for the stochastic bath used to keep the 
system at stationarity. In particular our thermostat is able to equilibrate the system 
also when collisions are elastic [19]. This happens because, in addition to a random 
driving, our thermostat acts on the particles also through a finite drag, in such a way 
that the temperature of the thermostat Tf,, different form the kinetic temperature of the 
fluid Tg < Tfj, is always well deflned. 

A remarkable feature of our model, related to the kind of thermostat we use, is the 
finite extent of velocity correlations in space. Indeed, the characteristic shape we find 
for transverse and longitudinal velocity structure factors allows us to define two non- 
equilibrium correlation lengths, ^ and ^i, which are known functions of the kinematic and 
longitudinal viscosities, respectively. This means that, instead of sampling trajectories 
of the system, out of equilibrium an average of static observables is enough to measure 
transport coefficients. The fiattening of the velocity structure factors at equilibrium 
clearly results in our formulas from the vanishing of the velocity correlations amplitude, 
which is proportional to Tb — Tg. In that case, access to transport coefficients is only 
possible through the study of the dynamical structure factors. 
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Finally, we study the behaviour of such coherence lengths at different packing 
fractions. We observe a significantly growth of the relative extent of correlations ^/Xq 
and ^i/Ao with the packing fraction, where Aq denotes the mean free path. This is 
in agreement with the slowing down of the dynamics also observed in dense granular 
fluids [20]. 

The paper is organized as follows. In Section [2] we discuss the model, the 
corresponding hydrodynamic equations and the specific forms of noise used. The k 
dependence of the eigenvalues of the linearized hydrodynamic matrix is also studied in 
order to check the stability and bounds of the linear approximation. In Section [3] we 
present a complete study of the static and dynamical structure factors for the velocity 
and density modes, comparing analytical predictions and numerical results. In section H] 
a discussion about the main findings is presented, with conclusions and perspectives for 



future work. In Appendix A and Appendix B are reported, respectively, formulas for 



the transport coefficients and some details on the noise terms. 

2. Microscopic model and fluctuating hydrodynamics 

We consider a system of N inelastic hard spheres in d dimensions with mass m and 
diameter a. Particles are contained in a volume V = L'^, with L the linear size 
of the system. We denote by n = N/V the number density and by the occupied 
volume fraction, (in two dimensions = N7r{a/2y/V). The particles undergo binary 
instantaneous inelastic collisions when coming at contact, with the following rule 

v. = v:-ii±^[(v:-v;).^]^ (1) 

where Vj (vj) and v^ (v') are the post and pre-collisional velocities of particle i (particle 
j), respectively; a G [0, 1] is the restitution coefficient (in the elastic case a = 1), and 
<T is the unit vector joining the centers of the colliding particles. 

In order to maintain a stationary fiuidized state, an external energy source is 
coupled to every particle in the form of a thermal bath [12] . In particular, the motion 
of a particle i with velocity Vj is described by the following stochastic equation 

mv,(t) = -7bV,(t) + ^,^,(t)+F,. (2) 

Here 7^ is a drag coefficient (which defines the characteristic interaction time with 
the external bath, t^^ = 7^/772), ^^^(t) is a white noise with (^^^(t)) = and 
{^b,ia{t)ib,ji3{t')) = '2Th'jf,6ij6ai3S(t — t') (Greek indexes denote Cartesian coordinates), 
while Fj represents the action of particle-particle inelastic collisions. The effect of 
the external energy source balances the energy lost in the collisions so that an out- of - 
equilibrium stationary state is attained [T^]. In particular, let us stress that such energy 
injection mechanism acts homogeneously across the whole system, differently from other 
mechanisms where the energy is directly supplied only to a part of the system, as for 
instance for fluids under shear, or for systems in conctact with vibrating walls. 

The stationary state is characterized by two time scales and two energy scales: the 
time scales are Tb and the mean free time between collisions Tc] the energy scales are the 
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temperature of the thermostat Tf, and the granular temperature Tg = '"^™ ' < Th 
(equal sign holds only if a = 1 or if rf, ^ Tc). In particular, in the dilute limit, 
the granular temperature Tg satisfies the following equation in the non-equilibrium 
stationary state ^T 



T, = n- ^ /(^wi-°') y. (3) 



where A^ = ^JmlT\2'^~^dl(y and x(0) is the pair correlation function at contact. The 
model has a well-defined elastic limit a — ?■ 1, where the fluid equilibrates to the bath 
temperature, Tg = Tb. The viscous drag term — 7feVj in Eq. ([2]) models the interaction 
between each particle and the thermostat. It is important to observe that 7^ is not 
related to the transport coefficients of the granular fluid and is fixed as a model 
parameter. As mentioned above, 7;, introduces a time scale Tb in the system that rules the 
tendency of particle to relax toward equilibrium at temperature Tb. The characteristic 
time of collisions, Tc, in all our simulations will be kept much smaller than Tb'. for this 
reason Tc is considered the microscopic time-scale of our system since it dictates the 
smallest scale of relaxation toward the non-equilibrium stationary state. In particular, 
in the coarse-grained hydrodynamic description to be introduced below, we will take 
care of comparing the characteristic decay time of different hydrodynamics modes with 
Tc, to verify the presence of a sufficient separation of scales. 

In the following we will present a thorough numerical analysis of model (jj]), using an 
event-driven molecular dynamics algorithm [21]. In particular, we will consider periodic 
boundary conditions in d = 2 dimensions. The fixed parameters of the simulations are 
m = 1, a = 0.01, Tb = 1 and 7^ = 1. The packing fraction is varied by changing the 
seize of the box, and we consider systems with G [0.1,0.5]. The simulation data on 
static structure factors are obtained for a system of A^ = 10000 particles, averaged over 
about 100 realizations, whereas the results on dynamical correlators are obtained with 
samples of iV = 1000 particles, averaged over about 4000 realizations. 

2.1. Linearized hydrodynamics 

Because we are interested in the behavior of large-scale spatial correlations in our system, 
we introduce here the coarse-grained hydrodynamic fields n(r, t), u(r, t) and T{r,t) as 
follows: 

n(r,t) = ^5(r-r,(t)), 

i 

n{r,t)=^J2^MS('^-r,{t)), (4) 

i 

_, , 2m -^-^ vf (t) ^ , , ,, 
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The hydrodynamic equations for the fields (jl]) can be derived for the model ([2]) following 
a standard recipe [221 



dtn{r, t) 


-V-(n(r,t)u(r,t)) 






dtu{r, t) + u • Vu = 


--V-n-^u(r,t) 
p m 




(5) 


9tT(r,t) + u- Vr = 


-—(v-j + n-.vu)- 

nd 


m 


-T(r,t)). 



In the above equations J and 11 are respectively the heat fiux and the pressure tensor, 
see details in Appendix B and 70 = (1 — a^)/2d. In the velocity equation the viscous 



drag term —'jbu/m has been inserted, while in the temperature equation three terms 
have been added: the sink term —T = —2'joUJcT{r,t) [25], where coc ~ A/T(r, t) is the 
collision frequency, takes into account the energy dissipated by inelastic collisions, while 
the terms 2'yb{Tb — T)/m represent the energy exchanged with the thermostat. 

Eqs. ([5]) give a fair description of the mesoscopic degrees of freedom of a granular 
fiuid as long as a proper separation of space and time scales is verified between those 
degrees of freedom and all the microscopic ones which are projected out. This condition 
is, of course, not always satisfied [261 [23 |28], but is not prevented in principle and is, 
indeed, realized in many experiments or simulations [291 EOl EH EQl 132] . 

Eqs. (|5|) can be linearized around the stationary homogeneous state, where 
the hydrodynamic fields take the values n = n, T = T and u = 
0. A system of linear differential equations for the fluctuations 5a(k, t) = 
{5n(k, t), 5T(k, t),M||(k, t),Mj_(k, t)}, with Sa = a— a, can be considered, with the Fourier 
transform deflned as 

5a(k, t)= f dr Sa{r, t)e-'^-', (6) 

and with M^(k, t) and M||(k, t) respectively the shear and longitudinal modes, namely 

U||(k) = A;-u(k) 

u^{k) = h-n{k), (7) 

k± being a unitary vector such that k± ■ k = 0. The system in Eq. ([5]) in Fourier space 
becomes 

Sa(k,t) = M{k)5a{k,t), (8) 

with the dynamical matrix 

/ tkn \ 

joUcg{n)Tg/n SjqUc + Dxk'^ + 27;,/m i2kp/dn 

ikv'^/n ik'p/ pTg Pik"^ + 7b/m 

\ Q Q vk"^ + '-fb/m J 



m{k) 



(9) 



where p = nm, Dt = 2i%/nd is the thermal diffusion coefficient [k is the heat 
conductivity), while v and vi are the kinematic and longitudinal viscosity respectively. 



Formulas for all parameters and transport coefficients are given in Appendix A There 
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Figure 1. Eigenvalues spectrum X{k) of the dynamical matrix (jH]) calculated for 
(f) — 0.5 and a = 0.6, namely in a strongly inelastic regime. The eigenvalues 
arc normalized on the Enskog collision frequency lJc so that on y axis the ratio 
Tc/T{k) = X{k)/ujc between the microscopic time-scale of the model and the decay 
time of each mode can be directly read. Inset: zoom of the spectrum at small k. 



we refer to the Enskog theory for dense elastic hard spheres (EHS) [33], which provides a 
good approximation, as observed in [18l[34]. The following sections are devoted to show 
how the viscosities u and ui can be obtained as fit parameters of static and dynamical 
correlations. Such results will be compared with the dense EHS predictions, finding 
good agreement. 



2.2. Spectrum of the hydrodynamic matrix and separation of time-scales 

We analyze the eigenvalues of M(/c) in order to study the linear stability of the model and 
to characterize the range of validity of time-scales separation required by hydrodynamics. 
From the expression in Eq. ([9]) we learn that the shear modes are decoupled from all the 
others, and the typical time-scale of their decay simply reads as r_L = (7fe/m + z//c^)~^. 
To obtain the typical time-scales for the fluctuations of the other hydrodynamic fields 
we solve, numerically, the equation Det(M(A;) — A(/c)I) = 0, where I is the identity 
matrix. The eigenvalues spectrum thus found is shown in Fig. [1], for parameters of the 
model a = 0.6 and (p = 0.5. All parameters are calculated according to the formulas 



reported in Appendix A 



The study of hydrodynamic eigenvalues shows that, even in the case of a quite 
inelastic and dense regime, a range of scales where mesoscopic relaxation times are larger 
than the microscopic times exists and a hydrodynamic description can be attempted. In 
particular, comparing the eigenvalues \i{k) with the Enskog frequency Uc (which will be 
verified, below, to be a good estimate of the real frequency in simulations), see Fig. [H 
we find, for each eigenvalue, the interval of values of k where Xi/uJc < 1 is fulfilled. 
The subscript i indicates eigenvalues related to heat {i = H), sound {i = ±) and shear 
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{i =-L) modes. The range is narrower for the case of the mode dominated by temperature 
fluctuations (denominated "heat mode" in the Figure, red thin curve), while it is larger 
for shear modes (green dot-dashed lines) and modes dominated by longitudinal velocities 
and density fluctuations (here referred to as "sound modes" , black and blue thick curves, 
continuous for real part and dashed for imaginary part). 

A general observation is that eigenvalues never have a positive real part, i.e. no 
instabilities are found, thanks to the presence of the external bath. Moreover, for the 
whole range of k studied, the two "sound modes" - as usual - are complex conjugate, 
i.e. they propagate with a /c-dependent velocity. A negative real part is always present, 
determining overall damping. At small values of k one can identify a sound velocity c by 
observing a linear relation Im[A-|-] ~ ck; instead the dispersion relation becomes strongly 
nonlinear at large k. Surprisingly, the study of sound eigenvalues shows that they have 
bifurcations at very small k < ki and very large k > k2 wave-numbers, becoming in both 
limits pure real numbers, i.e. losing their propagating behavior. Those bifurcations are 
due to the external damping, ruled by 7^: indeed in the hmit 76 — )■ (keeping finite 
jbTb), they disappear and the spectrum studied in [HI [20] is retrieved. Moreover, we 
find that the wavevector ki where the first bifurcation occurs, moves towards smaller 
values of k when the dissipation is increased, namely when a is decreased (at fixed 0) or 
when is increased (at fixed a). Looking closely to the bifurcations, we see that one of 
the eigenvalues for "sound modes" approaches for A; — )■ (representing total number 
conservation) and the other tends to — 7fe/m (as the shear one), see inset of Fig. [H The 
bifurcation at large k is perhaps non-physical, as it always falls out of the hydrodynamic 
range. For the eigenvalue of the heat mode we find that Xnik = 0) = —{3'joUc + 2jh/m). 
In the numerical setup used below, we have kmin = ^n/L > ki, so that the purely 
exponential decay of sound modes with large waves is never observed. 

2. 3. Stochastic description with fluctuating hydrodynamics 

In order to fully account for the spatial fluctuations of the hydrodynamic fields and for 
the decay in time of such fluctuations we must add some noise terms to the linearized 
hydrodynamic equations: the basic assumption under fluctuating hydrodynamics is the 
same as for average (deterministic) hydrodynamics, i.e. a good separation of scales 
between hydrodynamic fields and microscopic degrees of freedom. In the linearized 
hydrodynamic equations the small scale fluctuations have been projected out, but their 
feedback on large scale fluctuations can be recovered by a proper addition of noise terms 
to dynamical equations: 

5a(k,t) = M(A;)(5a(k,t)+f(k,t). (10) 

A derivation from first principles of the noise f (k, t) is beyond our scope. A kinetic 
theory with fluctuations has been recently proposed in [35l [36], for the homogeneous 
cooling regime, which is very different from our case. A similar treatment has been 
realized, only for the shear mode, in a driven case (random kicks without damping, i.e. 
lb = 0) [37], showing that noise can be safely assumed to be white, at difference with 
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the cooling regime. In such a case, with an additional but reasonable assumption on the 
two-particle velocity autocorrelation functions (namely, that such functions have only 
components in the hydrodynamic subspace), it is found that the fluctuation-dissipation 
relation for the internal part of the noise is satisfied, as already assumed in [25l [T8] . 
Following those previous studies, we will consider valid such an assumption. Notice 
that it does not imply that fluctuation-dissipation relations will be satisfied by the whole 
noise, which is composed by internal as well as external contributions. In summary we 
write 

/ \ 

e^^{'k,t) + 2ik/nde'''{k,t) 

er(k,t)+2A;/per(k,t) 
\ er(k,t)+^fc/per(k,t)} / 

where the two sources of noise for the hydrodynamic fields fluctuations are put in 
evidence: the first is the external contribution coming directly from the thermal bath, 
namely the stochastic force ^^ of Eq. ([2j); the second is internal and enters through the 
constitutive equations for the heat flux and the pressure tensor. A detailed discussion on 



f(k,t) 



:iii 



noises is presented in Appendix B The external and internal noises are Gaussian with 
zero average. The variances of external noises can be obtained directly from Eqs. ([2]) 
and (H 



\k,t)e'%k',t')) 



^Ml5(t-t05(k + k') 

dn m 



(er(k,t)er(k',o) = (er(k,t)er(k',t')) 



1 2^bTb 



6{t-t')6{k + k'),{12) 



n m 
while the variances of the internal contributions are obtained by imposing the 



fluctuation-dissipation theorem (see Appendix B for details): 



(^'"(k, t)^*"(k', t')) = 2KT^6{t - t')6{k + k') 
(^r(k, t)er(k', f)) = 2nmuiT,6{t - t')5{k + k') 
(er(k, t)er(k', t')) = 2nmuTg6it - t')Sik + k'). 
The internal and external noises are uncorrelated 



(13) 



(^^(k, t)^"(k', t')) = (er (k, t)c(k', t')) = (^r(k, t)^r(k', to) = o. (i4) 

The hydrodynamic analysis of model ([2]) consists, then, in solving the system of 
coupled linear Langevin equations fITOl) . In particular we are interested in finding the 
explicit forms of the static and dynamical structure factors, respectively 



Sab{k)= lim-(5a(k,t)(56(-k,t)). 



and 



SahiKu) 



Sab(k,t)e 



-iu)t 



(15) 



(16) 
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where 






0.1- 




Figure 2. C±{k,t) for a = 0.6 and 4> = 0.5 and several values of ka measured 
in numerical simulations, together with exponential fits (continuous curves). In the 
inset the quantity ^{k) obtained from the fits is reported, rescaled with the collision 
frequency uJc, together with the parabolic fit via the formula j{k) = 7b + vk^ 
(continuous red line). 



Sab{'k,t)= lim -{6a{k,t' + t)Sb{-k,t')). 



(17) 



3. Out-of-equilibrium correlations: Static and Dynamical structure factors 

It is well known that spatially extended correlations develop in the non-equilibrium 
stationary state of a driven granular fluid [HI [20]. In particular, the velocity correlator 

(u(k) ■ u(— k)), where the average (. . .) is taken over noises, can be written as 

(u(k) ■ u(-k)) = (M||(k)M||(-k)) + (Mx(k)Mx(-k)), (18) 

with the cross terms (n_|_(k)n||(— k)) = 0. The two terms on the right of Eq. flTSj) can 
be studied separately. 



3.1. Shear modes 

According to the matrix ([9|) the shear modes are decoupled from the others in the linear 
approximation and their dynamics obeys a simple Langevin equation 

Mx(k, t) = -(76 + iyk^)u±ik, t) + er(k, t) + ik/n ei"(k, t), (19) 

where, from now on, we put ?7i = 1 in all formulas. The effect of internal and external 
noises is equal to a single complex noise ^_L(k, t) = ^^.^(k, t) +ik/n C±{k., t) with variance: 



V-'{Uk,t)U{-Kt')) = -{T,^, + uTge)6{t-t'). 

n 



(20) 
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From this, the equal-time correlator {u±(\<.)uj_{—k.)) can be easily calculated: 

nSAt) ^ N-H\n4>^)l^) ^ ^-^t^ ^ T, . fl^, (21) 

with ^^ = 1^/76. The meaning of the above equation is clear: the small length-scale 
physics, which depends on the inelastic collisions, is related to the granular temperature 
Tg, while at large distances there are correlations with finite amplitude T^ — Tg and 
extent ^. Notice that in the limit 7^ — )■ 0, keeping 7fcT;, finite, from Eq. fl2Tl) one obtains 
the result of Ref. [TS], where S±{k) ~ l//c^ and long range correlations are observed. A 
similar power law behaviour is also observed in molecular fiuids under shear [7]. The 
different result we obtain in our model is due to the intrinsic cut-off introduced by the 
viscous drag 7b > 0. Indeed, the finite extent of correlations is even more clear when 
Eq. fl2T]) is written in real space, yielding the spatial correlation function G±{r), which 
reads: 

nGx(r) = Tg6^'\T) + (T, - Tg)^-^^, (22) 

where Ko{x) is the 2nd kind modified Bessel function that, for large distances, decays 
exponentially 

'<«(r/0 « y|(^. (23) 

At equilibrium, namely when collisions are elastic and Tg = Tb, equipartition between 
modes is perfectly fulfilled and the structure factor becomes fiat, i.e. S±{k) = T^. 
Differently, in the granular case, where Tg ^ Th, equipartition breaks down and from 
Eq. fl^ we have that S± — )■ T^ for small k and S± — )■ Tg for large k. We see here that out 
of equilibrium the quantity ^ = \/i^/% measures the range of static correlations of the 
vorticity field. The behaviour described above is in good agreement with experimental 
results obtained for driven granular fiuids, as reported in [38] and, more recently, in [32] . 
From Eq. (fT9l) we also find that fiuctuations decay exponentially 

{u±{k, t)u^{-k, 0)) ~ S^{k) e-(^^+'^'=')*, (24) 

with a characteristic time r(/c) = 7(/i;)~^ = (7^+///^^)"^. Such a behavior is also observed 
for elastic fiuids, with the only difference that in that case S±{k) = T^/n is constant. 
The length-scale ^ = \/i^/jb can be therefore always connected to dynamical properties 
of the system. What is peculiar of the out-of-equilibrium regime is that the so defined 
^ also represents the extent of correlations of the vorticity field, thus establishing a 
remarkable link between static correlations and dynamical ones. 

In Fig. [2] are reported the correlators C±{k,t) = {u±{k,t)u±{—k,0)) / {\u±{k)\'^) for 
different values of k and the same packing fraction measured in numerical simulations. 
Let us notice that the decay of {u±{k,t)uj_{—k,0)) is always exponential: we can 
therefore a-posteriori support the validity of linear hydrodynamics, from which the 
Langevin equations for shear modes is obtained. By interpolating with a parabola the 
characteristic time T{k) as a function of k, see inset of Fig. ([2]), we obtain the shear 
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Figure 3. Shear modes structure factor nS±{k) at different packing fractions for 
a = 0.6. Full curves arc drawn by inserting into Eq. (|2ip the values of v obtained from 
the decay of dynamical correlators. 






Statics 


Dynamics 


dense EHS 


0.1 


0.0062 


0.0072 


0.0111 


0.2 


0.0036 


0.0041 


0.0052 


0.3 


0.0025 


0.0031 


0.0039 


0.4 


0.0026 


0.0027 


0.0037 


0.5 


0.0020 


0.0028 


0.0041 



Table 1. Comparison of the shear viscosity obtained from statics (fit via Eq. (PTjl ). 
dynamics (fit via Eq. (|24p ) and dense EHS approximation. 



viscosity u. Let us stress the deep connection between statics and dynamics in the 
out-of-equihbrium regime: inserting the values of u obtained from the dynamics into 
Eq. ( 12T|) . we find curves that well superimpose the numerical data for {u±{k)u±{—k)), 
(see Fig. [3]). The values of u obtained from the dynamics of shear modes can be 
independently obtained as fit parameters of S±{k) via Eq. ( 12T]) . The values of z/ obtained 
with the two different procedures are compatible, as can be seen from Tab. [1] They are 
also reasonably close to the dense EHS predictions, presented in the same Table. 



3.2. Longitudinal modes 

3.2.1. Static correlations The same considerations discussed above for shear modes 
also hold for the other hydrodynamic modes, which are coupled each other. In order to 
study their behaviour we have to take into account all the elements of the dynamical 
matrix. In particular, the matrix of static structure factors S(k) with elements Satik), 
is obtained solving the following linear system: 

M{k)S{k) + S\k)M\k) + C{k) = 0, (25) 
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where the matrix of noises C{k) is such that: 

V^' (f (k, t) ® f (-k, t')) = C(k)5(t - t'). (26) 



Here 

C(k) = diag 



2T 2 2 

0, -^ (4Tfe7fe + 2DTTgk'') , - (Tb7, + uiT^k'') , - [T,^, + uTgk'') 
dn ^ ' n^ ^ n ^ ' 



(27) 



where diag[x, y, z, vo\ denotes a diagonal matrix with elements x, y, z, w. 
The expression of the longitudinal structure factor, 

nS^\{k) = N~'{\ui{k)\^), (28) 

turns out to be the ratio between two even polynomial functions of the 6— th order in 
k: 

11^^ ~ S', + S',k^ + S',k^ + S',k^- ^ ' 

In the above expression eight constants have been introduced, which depend in a 
complicate manner by all the parameters of the system. Let us focus here on the 
asymptotic behavior of S*]] (k) at large and small values of k. From Eq. ([3]) there follows 
the relation djoUJc/2 = jb{Tb — Tg)/Tg, which allows us to recast the series expansion 
around fc = of the expression in Eq. (129|) in the form 

n5||(fc ^ 0) ~ T, - (n - Tg)^fe + 0{e), (30) 



with 



.2 _ ^r _ 1 

Si 



^^ ^ Ibll / Jp^ ^ 2g{n)p 

raTg(7ft + 7oa;c)(27f, + 37oa;c) \d^nTg 3d 



(31) 



lb lb 
Up to 0{k'^) the expression in Eq. fl30l) is equivalent to 

^^llW = r. + 3^^> (32) 

namely a form analogous to the structure factor we found for the shear modes. In Fig. H] 
we show the static structure factors for different packing fractions. Again, notice that in 
the limit 7(, — )■ with ibTh finite, the behaviour S\\{k) ~ 1/A;^ found in [18] is recovered. 
Albeit the above expression is in principle only valid for low k values, it captures also 
the large k limit, when fine oscillations are disregarded. Indeed, expanding Eq. ( !29ll for 
large k values we find: 

nS\\{k^oo)^Tg+ ^^'~^^'\ (33) 

Such a discussion shows that even for longitudinal modes the viscosity vi is related to a 
finite correlation length, measurable from static velocity correlations when the system is 
out of equilibrium with Tg < Th. The behavior of that length when the packing fraction 
is increased will be discussed in the last section. 
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Figure 4. Longitudinal modes structure factor nS'n (fc) at different packing fractions, 
for a = 0.6. Full curves are drawn inserting into Eq. (j32p the values of ^i, defined in 
Eq. (PT|) . obtained from best fits of the dynamical structure factors. 



3.2.2. Dynamical correlations Dynamical correlations for longitudinal modes are less 
simple than those for the shear mode, since they are given by a superposition of different 
(real and imaginary) exponentials. The dynamical structure factors are obtained by 
solving the equation of motion ( ITOl) . which, in the frequency domain, reads as 

M(k, w)5a(k, uj) = f (k, w), (34) 

where 

M(k,uj)=zujl-M{k), (35) 

with I the identity matrix, 

/oo 
dt 5a(k,t)e-''^\ (36) 

-oo 

and 

V~\{{k,uj)^i{-k,u')) = C{k)6{uj + uj'), (37) 

C{k) being defined in Eq. (E?]). Muhiplying Eq. §S) on the left by M-\k,uj) and on 
the right by 6si^{—k, —u) (where X'^ denotes the transpose of X) and averaging over 
the noise, we obtain the matrix of dynamical structure factors 

S{k,uj) = V~\M~\k,u)f{k)6si^{-k,-u)) = M-\k,u)C{k)[M^{-k,-u)]-\ (38) 

where in the last equality we have used the Hermitian conjugate of Eq. (IMl) and the 
relation (13 
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The dynamical structure factors Snn{k,u) and Sii{k,u) take the exphcit forms 



2l2 



Snn{k,uj) = n k 



where 



/ [u' + (27. + Dt^ + 37oc^c)^] (^ + ^; 



\nTg ) \ nd """ nd ) 



\ 



detMP 



(39) 



detM|2 



,2/o,^,, , n l2 , .,,2 , o,^A , l2//o,^ , o,^,,N.,2 ^WPTO^c , o 



-uj\?>-1qUJc + I^tA; + vk' + 37fe) + ^ (27fe + 37oWc)«r " 



+ < cj" - a; 



and 



7b(27fe + 37oa;c) + k' Driiyk' + 7.) + 1^(27. + 37oWc) + 



n 



2p2 



+ v^Drk' 



S\\{k,uj) 



UJ 



v?k'^ 



*^nn[k, OJ)- 



dn^mTg '^ 



(41) 



1 2 



(40) 
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Sim 


dense EHS 
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160 


181 


181 


= 0.3 




a = 0.6 


a = 0.8 




Sim 


dense EHS 


Sim 


dense EHS 


T 


0.093 


0.0829 


0.125 


0.1185 


UJc 


79 


83 


85 


85 


= 0.1 




a = 0.6 


a = 0.8 




Sim 


dense EHS 


Sim 


dense EHS 


T 


0.212 


0.2055 


0.286 


0.2820 


Wc 


26 


25 


29 


28 



Table 2. Comparison of theoretical predictions of Eqs. ([3]) and (jA.2[) and numerical 
results for Tg and lOc- 

All the coefficients appearing in these expressions can be evaluated within the dense 
EHS approximation, so that Eq. ( 14T|) can be used in order to obtain z// and Dt from the 
numerical data. Clearly the writing of physical quantities, such as v^ = {dp/dn)T or p 
itself, with dense EHS formulas leads to systematic errors. The amplitude of this error 
can be estimated for instance by comparing the theoretical prediction of dense EHS 
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Figure 5. Dynamical structure factor nS^^{k,uj) at fixed paclcing fractions (j) and 
different momenta. Identifying a characteristic time Tp^ak with the frequency of the 
maximum of the curve we see tliat Tc/Tpgak = w/ojc grows with the wave- vector, namely 
longitudinal modes of higher momentum decay faster. Continuous lines show the best 
fit results via Eq. (|4T|) . with Dt and i^i as fitting parameters and all other coefficients 
fixed with the dense EHS predictions. 



collision frequency with the collision frequency measured in simulations. In Table |2] are 
listed both dense EHS and numerical collision frequencies at different packing fractions. 
Again it is found that the dense EHS prediction is quite good. 

In order to obtain the longitudinal viscosity and the thermal diffusion coefficient, we 
fit our numerical data for the longitudinal modes using Eq. ( 1¥T|) . where all parameters 
but ui and Dt are fixed to the dense EHS values and Tg is the one measured in 
simulations. In Fig. |5] S\\{k,u) is shown for different values of k and fixed = 0.5, 
together with the best fit curves. The values of ui and Dt so obtained, together with 
those computed within the dense EHS approximation, for different values of </> and a, are 
reported in Table [3] and within errors are found independent of k. This fact represents 
an a-posteriori check that we are in the regime of validity of linearized hydrodynamics. 
Indeed, in Eqs. (HTj) and (140|) z// appears as fc-independent variable. Only at low packing 
fraction we observe a dependence on k. This is perhaps due to diluteness, which implies 
too large mean free path or mean free time with respect to mesoscopic scales. 

In Fig. |6] we also show the time decay of the dynamical structure factor. It can 
be appreciated, in the time domain, the superposition of different real and imaginary 
exponentials, which determines a mix of damping and propagation. 

The dynamical structure factor at fixed k and different packing fractions is reported 
in Fig. [3 it is remarkable the observation of a time-scale, individuated by the peak 
frequency of S\\{k,u), which increases as the packing fraction is increased. Such 
behaviour is consistent with the observation, discussed in details below, of a growth 
of the correlation lengths defined above, together with the packing fraction. 
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Fit Results: ka = 0.2 


0.039 
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ka = 0.3 
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0.019 


0.018 


ka = 0.4 


0.018 


0.013 


0.0096 


0.016 


dense EHS 


0.048 


0.012 


0.056 


0.014 



Table 3. Comparison of theoretical predictions of Eqs. (|A.6IIA.10| and fit results via 
Eqs. dH]) for Dt and i^i. 



We conclude this section by stressing the remarkable agreement between numerical 
S\\{k) data and the expression of Eq. (132|) . see Fig. HI with ui measured from dynamics 
and the renormalization term entering the definition of z/^* (see Eq. (l3T]) ) calculated 
within the dense EHS approximation and reported in Table [3l 



4. Summary and conclusions: transport coefficients and non-equilibrium 
correlation lengths 

In conclusion, we have studied both static and dynamical correlations for hydrodynamic 
fluctuations of the velocity and density fields: indeed we recall Eq. fHT]) which gives 
a direct relation between density and longitudinal velocity structure factors. A 
main comment concerns the good success of analytical predictions: comparison with 
simulations shows a fair agreement up to = 0.5, with values for most of the parameters 
directly given in the dense EHS approximation. This signals a success of the three 
main ingredients: 1) dense EHS theory for transport coefficients, working even at quite 
high densities, 2) scale separation and granular hydrodynamics, and 3) prescription 
of Eqs. ( IT2IIT4I) for the hydrodynamic noise. The last ingredient is perhaps the most 
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Figure 6. Dynamical structure factor nS^^{k,t) at fixed packing fraction (p and 
different momenta. The observed oscillations are in agreement with the eigenvalues 
spectrum in the interval of momenta considered: the eigenvalues of sound modes are 
complex conjugate, thus producing an oscillatory relaxation. 
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Figure 7. Dynamical structure factor nS\i{k,uj) at fixed momentum k and different 
packing fractions. Identifying a characteristic time Tpeak with the frequency of the 
maximum of the curve we see that Tc/Tpeak = w/wc grows when the packing fraction 
is lowered, namely sound modes decays faster compared to the microscopic time-scale. 
Continuous lines show the best fit results via Eq. (HI]), with Dt and i^i as fitting 
parameters and all other coefficients fixed with the dense EHS predictions. 
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interesting, if one considers how hard is the characterization in simple terms of non- 
equihbrium systems. Moreover, the assumption of white internal noise is less obvious 
than in the case of simple random kicks without drag [37]: in that case, the absence of 
drag (7b = 0) implies a relaxation time at large scales, fc — )■ 0, which diverges, making 
a solid base for assuming fast the relaxation of noise due to microscopic degrees of 
freedom. In our model, in principle, the drag could be large enough to make even large 
scales fast, making difficult to define hydrodynamic fields. Recent experiments show 
that, even if that could be a realistic situation, this does not change dramatically the 
qualitative behavior of structure factors [32] . 

One of the main results of our study is the presence of spatial order in the form 
of non-equilibrium velocity correlations. This should be related to the slowing-down of 
the dynamics with increasing packing fraction in granular systems [39l [20], and with 
the existence of a time-scale growing with the density, as observed in [TT]. As discussed 
in Sec. [3l in our model the competition between different relaxation mechanisms given 
by the kinematic and longitudinal viscosities u and ui and the thermostat damping 7^ 
give place to a couple of length-scales characterizing non-equilibrium structure factors: 
^ = ^/^/ib and ^i = \fi^fj^h- The first trivial observation is that allowing 7;, — )■ 
such lengths diverge. It means that, according to the prediction of [18], the largest- 
scale correlations are always equal to the size of the system, which is the largest size 
available. In our model there is a cut-off on such correlations imposed by the viscous 
drag 7fe > due to the interaction with the thermal bath. The existence of such a 
cut-off, which represents a fixed parameter at different packing fractions, allows us to 
show that by increasing the packing fraction the extent of correlations is eflectively 
increased. While the absolute value of ^ is only slightly enhanced when is increased, 
as also observed in experiments in [38], we stress that, in order to appreciate the physical 
meaning of such length at different densities, one has to compare it with the microscopic 
relevant spatial scale in the system, which is given by the mean free path of the particles 
^0 = "\Ai"cr/(80x(0)) cind which also changes with the packing fraction. Then one 
finds that .^/Aq is remarkably increased at high densities, as can be seen in Fig. [8], and 
also in recent experiments [32]. We may summarize the observed phenomenon saying 
that the higher is the packing fraction the higher must be the number of intermediate 
scattering events between two different particles in order to decorrelate their velocities. 
This scenario could be reflected in transitions of dynamical origin, at higher packing 
fractions [iQl SH iSl [39] . 
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Figure 8. Correlation lengths ^ and £,i rescaled with the mean free path Aq, for 
different packing fractions and several values of the restitution coefficient. The extent 
of correlations grows exponentially with the packing fraction. Notice also that, at fixed 
(j), ^ and ^i are larger for higher values of a. However the amplitude of correlations, 
ruled by Tb — Tg, is smaller, see Eqs. (f2T|) and (|30l). 



Appendix A. Dense EHS formulas for hydrodynamic coefficients 

The definitions of tlie parameters entering the matrix M are 



jj. 



UJr 



dp 
dp 






T„ 



g{<t>) = 21 + 



'9x(0) \ 

In the dense EHS approximation the pressure p can be written as 

VtdXncr'^ 1 



p{n) =nTg[l 



a 



2d 



which in d = 2 reads as 

/I 

p{4>) 



vro" 



,0Tg(l + 0x(0)(l + a)), 



(A.l) 

(A.2) 
(A3) 

(A.4) 
(A5) 



where has been made use of the relation between packing fraction and density n = 
40/(7rcr^), and of the definition Vtd = 27i'^^'^/T(d/2) with d = 2. Notice also that we have 
taken into account the correction due to the inelasticity, as given in [13] . There are dense 
EHS formulas also for the dependence of the diffusion coefficients on the packing fraction 
and the granular temperature. In the following such formulas are written for a 2d system. 
In this case, we use the Verlet-Levesque approximation for the pair correlation function 
at contact: x(0) = (1 - 70/16)/(l - 0)^. 
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(A.6) 



(A.7) 



(A. 



(A.9) 



(A.IO) 



Appendix B. Noises 

In this appendix we present a detailed discussion of the noise terms appearing in the 
fluctuating hydrodynamic equations ( ITOl) . Let us start from the external noises, which 
can be simply obtained from Eqs. ([2]) and (j4]). We find that the noise contributions to 
the equations for the velocity and temperature fields are, respectively 



r{r,t) = -Y,^,Amr-rm 



n 



2m 



0'^{r,t) = —J2Mt)-^,Amr-rm- 



These are Gaussian noises with variances 



(B.l) 
(B.2) 



(C(r,t)er(r',t')) 



ex,„/ ./^^ _ 1 '^^'''^'6^^6{t - t')6{r - r') 



{r,t)9^%r',t')) 



n m 
AmTg 2'ybTb 



S{t-t')Sir-r'). 



(B.3) 



dn m 

In order to describe the local spontaneous microscopic fluctuations of the fluid, we 
also include in our description internal conserved Gaussian noises 0*" and ^*", which 
enter the constitutive equations for J and 11, respectively 



kVT + e' 



n = pi — ?7 



Vu + fVu) 



-r/-C)lV-u + r 



(B.4) 
(B.5) 
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where k is the heat conductivity, p is the local pressure, 1 the unit tensor, rj the shear 
viscosity and ( the bulk viscosity. The amplitudes of such noises are obtained from the 
fluctuation-dissipation theorem |H O HI] 

(Q(r,t)e;'^(r',t')) = 2TMSa,5ps + S^sS^,) + U - ^v) S^pS^sW - t')Siv - r') 
{e^{v,t)ef{r\t')) = 2KTXAt - t')^{^ - r'). (B.6) 



Notice that, for granular systems, Eq. f lB.4p should be modified, adding the term — /iVn 
on the rhs, which takes into account the contribution to the heat current due to density 
gradients |15] . However, the transport coefficient /i is very small for driven systems |16] , 
and therefore we have neglected this contribution. 

Linearization 

In order to obtain the linear approximation of the hydrodynamic equations, we consider 
how a homogeneous fluctuation of the temperature relaxes around its stationary value 
Tg. Therefore we start by linearizing Eq. (jS]) around Tg, with ?T,(r, t) = n, u = 0, T(r, t) = 
T{t), and obtain 

5f{t) = -(27b/m + 3^oUc)5T{t). (B.7) 

Notice that the quantity —{2%/m + S'-y^Uc) coincides with the heat mode eigenvalue 
XHik = 0). 

Next, taking the linear terms around the non-equilibrium steady state in Eqs. (jS]), 
introducing the external noises fIB.ip and flB.2p . and using Eqs. f lB.4p and f IB.Sp . we 
obtain 



dtSn{r,t) = — nV-u(r,t) 

dtu{r, t) = - - Vp(r, t) + z/V^u(r, t) + (ui - z/) VV ■ u(r, t) 

- ^u(r, t)--V- r (r, t) + r{r, t) (B.8) 
m p 

d,6T{r, t) = ^V^5T(r, t) - ^ V . u(r, t) - l^^^^MllSn 
nd nd n 

- 2^6T{r, t) - 3-fou, 6T{r, t) - —y ■ 6>^"(r, t) + ^^(r, t), 

m nd 

with pu = 7] and pui = 2ri{d — l)/d + (. 

The velocity field u can be split in longitudinal and transverse components u = 
u/ + ux, where V ■ u^ = and V A uy = 0. With such a decomposition Eqs. (IB.SP can 
be written 



dt6n{r,t) = — ?7,V-u(r, t) 
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dtui{r,t) = - -Vp(r,t) + z/,VV ■ uKr,t) - ^u^r.t) - -V ■ C(r, t) + ^r(r, t) 
p m p 

atux(r, t) = uVW{r, t) - ^u^(r, t)--V- ^^(r, t) + ^^(r, t) (B.9) 

m p 

dt6T{r,t) = ^V^6T{r,t) - %V ■ u(r,t) - 2^5T(r,t) - 870^;^ ST{r,t) 
na na m 

'yo(^c9{n)Tg ^^ 2 



5n-— V-0*"(r,t) + r"(r,t), 



where 



n nd 



(eia/3(r, t)ei:^5(r', t')) = 2Tgr){S^^Sps + S^sSp^ - 2S^p5^s)S{t - t')S{r - r' 

2r]{d-l 



(C/5(r,t)a;.(r',0) =2Z 



C + 



d 



6^p6^sS{t-t')6{r-r'), 

(B.IO) 



and 



(er(r,t)er(r',t')) = (er(r,t)r(r',t')) = -^^6{t-t')6{r-r').{B.ll) 

n m 



Finally, taking the Fourier transform of Eqs. (IB. 91 IB. 101 IB. Ill) we obtain Eq. (ITOl) and 
foUowings. 
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